#include <stdio.h>
#include <stdlib.h>
#include <time.h>

float max_error(int n, float *A, float *x, float *b);
void multiply(int n, float *M, float *v, float *result);
void sum(int n, float *a, float *b, float *result);


int main(int argc, const char* argv[])
{
   struct timeval ti, tf;
   long int sec, usec;
   //Início da contagem de tempo
   gettimeofday(&ti, NULL);
   
   
   float *A, *b;
   float *B, *g;
   float *x, *x_p;
   int n;
   int J_ORDER, J_ROW_TEST, J_ITE_MAX;
   float J_ERROR;
   int i, j;
   int n_iter;
   float error;
   FILE *f;
   
   f = fopen(argv[1], "r");
   fscanf(f, "%d", &J_ORDER);
   n = J_ORDER;
   fscanf(f, "%d", &J_ROW_TEST);
   fscanf(f, "%f", &J_ERROR);
   fscanf(f, "%d", &J_ITE_MAX);
   
   A = (float*) malloc(n * n * sizeof(float*));
   B = (float*) malloc(n * n * sizeof(float*));
   
   b = (float*) malloc(n * sizeof(float));
   g = (float*) malloc(n * sizeof(float));

   x = (float*) calloc(n, sizeof(float));  
   x_p = (float*) malloc(n * sizeof(float));   
   
   for(i = 0; i < n; i++)
   {
      for(j = 0; j < n; j++)
         fscanf(f, "%f", &A[i * n + j]);
   }
   for(i = 0; i < n; i++)
      fscanf(f, "%f", &b[i]);
   
   for(i = 0; i < n; i++)
   {
      for(j = 0; j < n; j++)
         if(i != j)
            B[i * n + j] = -(A[i * n + j] / A[i * n + i]);
         else
            B[i * n + j] = 0;
      g[i] = b[i] / A[i * n + i];
   }      
      
   for(n_iter = 0, error = J_ERROR; n_iter < J_ITE_MAX && error >= J_ERROR; n_iter++)
   {
      for(i = 0; i < n; i++)
         x_p[i] = x[i];
      multiply(n, B, x_p, x);
      
      sum(n, x, g, x);
      error = max_error(n, A, x, b);
   }
   
   multiply(n, A, x, x_p);

   for(i = 0; i < n; i++)
      printf("x[%d] = %f\n", i, x[i]);
   printf("Iteracoes: %d\n", n_iter);
   printf("Linha %d: calculado = %f, esperado = %f.\n\n", J_ROW_TEST, x_p[J_ROW_TEST - 1], b[J_ROW_TEST - 1]);

   //Término da contagem de tempo
   gettimeofday(&tf, NULL);
   sec = tf.tv_sec - ti.tv_sec;
   if(tf.tv_usec < ti.tv_usec)
      sec--;
   usec = (tf.tv_usec - ti.tv_usec + 1000000) % 1000000;
   printf("Tempo decorrido: %ld.%06ld s.\n", sec, usec);
   return 0;
}

float max_error(int n, float *A, float *x, float *b)
{
   float *errors;
   float error;
   int i;
   
   error = 0;
   errors = (float*) malloc(n * sizeof(float));
   
   multiply(n, A, x, errors);
   for(i = 0; i < n; i++)
   {
      errors[i] = b[i] - errors[i];
      if(errors[i] < 0)
         errors[i] = -errors[i];
      if(error < errors[i])
         error = errors[i];
   }
   
   return error;      
}

void multiply(int n, float *M, float *v, float *result)
{
   float *res;
   float s;
   int i, j;
   
   res = (float*) calloc(n, sizeof(float));
   for(i = 0; i < n; i++)
      for(j = 0; j < n; j++)
         res[i] += M[i * n + j] * v[j];

   for(i = 0; i < n; i++)
      result[i] = res[i];
      
   free(res);  
}

void sum(int n, float *a, float *b, float *result)
{
   float *res;
   int i;
   
   res = (float*) calloc(n, sizeof(float));
   for(i = 0; i < n; i++)
      res[i] += a[i] + b[i];
   
   for(i = 0; i < n; i++)
      result[i] = res[i];
      
   free(res);  
}
